2 - Rasters in R
The first time you use a new version of R you will have to install all your packages. Try installing the raster and rgdal packages if you don’t have them already:
if(!require("raster")) install.packages("raster")## Loading required package: raster## Loading required package: spif(!require("rgdal")) install.packages("rgdal")## Loading required package: rgdal## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /home/mmann1123/.local/share/proj:/usr/share/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.if(!require("lwgeom")) install.packages("lwgeom")## Loading required package: lwgeom## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.10.2, PROJ 8.2.1if(!require("sf")) install.packages("sf")## Loading required package: sf## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUEif(!require("dplyr")) install.packages("dplyr")## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionif(!require("osmdata")) install.packages("osmdata")## Loading required package: osmdata## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyrightRaster data is stored in a format similar to that of a matrix of numbers, except you have added data about the coordinate references system (projection).
Let’s read in one raster and take a look
Examine Raster Data
Notice that the working directory doesn’t need to be a path directly below your files of interest:
CWD_apr = raster('./data/raster/CWD/cwd2010apr.tif')Get info about the raster object
CWD_apr## class : RasterLayer
## dimensions : 1120, 872, 976640 (nrow, ncol, ncell)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## source : cwd2010apr.tif
## names : cwd2010apr
## values : 0, 150.9694 (min, max)Plot the image:
plot(CWD_apr)
there are a variety of different plot options. Here are a few important ones
plot(CWD_apr,main = 'Title')
plot(CWD_apr,
col= heat.colors(10),
main='Water Deficit' ) # 10 designates the number of hex color codes that get generated
Now use ?plot to change col to a different preset color scheme
?plot## Help on topic 'plot' was found in the following packages:
##
## Package Library
## terra /home/mmann1123/R/x86_64-pc-linux-gnu-library/4.1
## raster /home/mmann1123/R/x86_64-pc-linux-gnu-library/4.1
## graphics /usr/lib/R/library
## base /usr/lib/R/library
## sf /home/mmann1123/R/x86_64-pc-linux-gnu-library/4.1
##
##
## Using the first match ...Raster Math
Rasters like matrices can make basic math functions very easy. Note that rasters can be used in basic operations like CWD_apr *2 or CWD_apr^2.
However teo rasters with different extents, cell sizes, or
projections CAN’T be used. You would need to reprojected and aligned
first, this can be done using the projectRaster() function,
run ?projectRaster for help.
CWD_apr_sq = CWD_apr^2
plot(CWD_apr_sq, main='CWD Squared')
Now lets add two different rasters. Use the raster() function to read in the CWD for May
Now let’s add two months of data values:
#CWD_may = raster('fill/in/the/path/file.tif')
CWD_A_M = CWD_apr + CWD_may
plot(CWD_A_M)
Let’s reproject that raster into a UTM zone 10 and write it back out
# Define the Proj.4 spatial reference
# https://spatialreference.org/ref/epsg/4326/
sr = "EPSG:4326"
# Project Raster
projected_raster = raster::projectRaster(CWD_apr, crs = sr)
# Write the RasterLayer to disk (See datatype documentation for other formats)
writeRaster(projected_raster, filename="./output/CWD_apr_UTM.tif", overwrite=TRUE)Raster Values
We can also directly set the cell values of rasters. For instance, it
is often useful to have a raster filled with a fixed value, such as
0. In this example we will use the setValues()
function to update the values of CWD_may to all be 0, and
store it in a variable called zeros.
zeros = setValues(CWD_may,0)
plot(zeros)
Note that zeros raster has the same extent, resolution,
and projection as CWD_may
print(zeros)## class : RasterLayer
## dimensions : 1120, 872, 976640 (nrow, ncol, ncell)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## source : memory
## names : cwd2010may
## values : 0, 0 (min, max)print(CWD_may)## class : RasterLayer
## dimensions : 1120, 872, 976640 (nrow, ncol, ncell)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## source : cwd2010may.tif
## names : cwd2010mayWe can also use boolean queries to update certain values. Here we
update certain values of CWD_may to 5000 where
CWD_may>150 is TRUE with the following
syntax:
CWD_may[CWD_may>150] = 5000
plot(CWD_may)
Note that we have now changed the values of CWD_may… if
we want to original values again, we will need to read them back in
again with
CWD_may = raster('./data/raster/CWD/cwd2010may.tif').
Challenge Question
In this section we will try to get the average water deficit for
California for the whole year of 2010. To do this we will need to
isolate the files we are interested in using - those ending with
.tif. We will use the function list.files()
which creates a list of all the files in a particular folder.
First lets get a list of all files in our raster folder:
list.files(path="./data/raster/CWD/")## [1] "cwd2010apr.tfw" "cwd2010apr.tif" "cwd2010apr.tif.aux.xml"
## [4] "cwd2010apr.tif.xml" "cwd2010aug.tfw" "cwd2010aug.tif"
## [7] "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml" "cwd2010dec.tfw"
## [10] "cwd2010dec.tif" "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"
## [13] "cwd2010feb.tfw" "cwd2010feb.tif" "cwd2010feb.tif.aux.xml"
## [16] "cwd2010feb.tif.xml" "cwd2010jan.tfw" "cwd2010jan.tif"
## [19] "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml" "cwd2010jul.tfw"
## [22] "cwd2010jul.tif" "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"
## [25] "cwd2010jun.tfw" "cwd2010jun.tif" "cwd2010jun.tif.aux.xml"
## [28] "cwd2010jun.tif.xml" "cwd2010mar.tfw" "cwd2010mar.tif"
## [31] "cwd2010may.tif" "cwd2010nov.tif" "cwd2010oct.tif"
## [34] "cwd2010sep.tif"Notice that .png and twf xml files are also listed. We can narrow this down by using a grep style pattern search.
Here we will only list files that have .tif in the
name:
list.files(path="./data/raster/CWD/",
pattern = '.tif')## [1] "cwd2010apr.tif" "cwd2010apr.tif.aux.xml" "cwd2010apr.tif.xml"
## [4] "cwd2010aug.tif" "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml"
## [7] "cwd2010dec.tif" "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"
## [10] "cwd2010feb.tif" "cwd2010feb.tif.aux.xml" "cwd2010feb.tif.xml"
## [13] "cwd2010jan.tif" "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml"
## [16] "cwd2010jul.tif" "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"
## [19] "cwd2010jun.tif" "cwd2010jun.tif.aux.xml" "cwd2010jun.tif.xml"
## [22] "cwd2010mar.tif" "cwd2010may.tif" "cwd2010nov.tif"
## [25] "cwd2010oct.tif" "cwd2010sep.tif"almost there… we still have problems with .tif.* files like tif.xml
or tif.aux.xml. We can limit these using what I will call an
ANTI-wildcard. “$”. '.tif\$' will tell the computer the
file needs to END in .tif. Let’s also store the path/file
name by setting full.names=TRUE
list.files(path="./data/raster/CWD/",
pattern = '.tif$',
full.names=TRUE)## [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
## [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
## [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
## [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
## [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"Let’s store that information in a variable called
CWDS:
CWDS = list.files(path="./data/raster/CWD/",
pattern = 'tif$',
full.names=TRUE)Challenge: In teams read in and sum all of these rasters, then divide by the total number of files.
# Hint
jan = raster("./data/raster/CWD//cwd2010jan.tif")
feb = raster("./data/raster/CWD//cwd2010feb.tif")
# etc
# sum them then divide by total number of files (need to include other months)
CWD_mn_2010 = (jan+feb ) / 2### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ### Finally, let’s write out our results:
# writeRaster(x=yourraster,filename="./output/CWD_mn_2010.tif")Raster Stacks
R also has the capacity to handle multi-band or multi-layer raster
‘stacks’. 
In this
case we will be using a stack to combine multiple files of the same type
(CWD) over time (months of 2010). As you will see there are multiple
advantages to handling data in stacks.
The only requirement is that all the images must have the same: - Image extent - Resolution - Projection
In this case our CWD data is already identical in shape and
projection. All we need to do is use a list of files like
CWDS that we created earlier.
CWDS## [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
## [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
## [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
## [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
## [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"and create a stack from the data:
cwd_stack = stack(CWDS)Now looking at the properties we can see that cwd_stack
has 12 layers.
cwd_stack## class : RasterStack
## dimensions : 1120, 872, 976640, 12 (nrow, ncol, ncell, nlayers)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## names : cwd2010apr, cwd2010aug, cwd2010dec, cwd2010feb, cwd2010jan, cwd2010jul, cwd2010jun, cwd2010mar, cwd2010may, cwd2010nov, cwd2010oct, cwd2010sep
## min values : 0, 0, 0, 0, 0, 0, 0, ?, ?, ?, ?, ?
## max values : 150.96939, 205.35623, 47.32375, 65.93313, 32.79562, 233.96126, 219.87500, ?, ?, ?, ?, ?This allows us to do operations much easier. For instance we can just run arithmetic functions on the entire stack:
CWD_mn_2010 = mean(cwd_stack)
plot(CWD_mn_2010)
We can also do ‘global’ statistics like finding the ‘mean’ of all cells for any given month.
cellStats(cwd_stack, stat='mean', na.rm=TRUE)## cwd2010apr cwd2010aug cwd2010dec cwd2010feb cwd2010jan cwd2010jul
## 45.5497388 157.9214946 1.6478612 5.2903873 0.1809682 176.2727036
## cwd2010jun cwd2010mar cwd2010may cwd2010nov cwd2010oct cwd2010sep
## 143.6983408 27.6639933 90.9696350 19.2020629 27.7737959 120.8165948Or other summary visualization like histograms on a month by month basis.
hist(cwd_stack)
Raster Vector Interactions
Raster cropping
Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case raster cropping and masking are useful for unifying the spatial extent of input data. Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.
We will use two objects to illustrate raster cropping:
- A SpatRaster object
srtmrepresenting elevation (meters above sea level) in south-western Utah - A vector (sf) object
zionrepresenting Zion National Park
Both target and cropping objects must have the same projection. The
following code chunk therefore not only reads the datasets from the
spDataLarge package. This package needs to be installed
directly from github using the following code:
# this has already been done for you, please don't run
#install.packages("spDataLarge", repos = "https://geocompr.r-universe.dev")We will also need to install/import the sf package to
handle vector data.
if(!require("spDataLarge")) install.packages("spDataLarge")## Loading required package: spDataLargeif(!require("sf")) install.packages("sf")
srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))Let’s take a look at the two files, note the use of
add=TRUE in the second plot command. This adds the new
object to the previous plot.
plot(srtm)
plot(zion, add=TRUE)
We use crop() from the raster package to crop the
srtm raster. The function reduces the rectangular extent of
the object passed to its first argument based on the extent of the
object passed to its second argument.
srtm_cropped = crop(srtm, zion)
plot(srtm_cropped)
plot(zion, add=TRUE)
Note that the extent of the strm raster has been cropped to
fit the zion polygon.
Related to crop() is the terra function
mask(), which sets values outside of the bounds of the
object are converted to NA, or a missing value. The
following command therefore masks every cell outside of the Zion
National Park boundaries:
srtm_masked = mask(srtm, zion)
plot(srtm_masked)
Importantly, we want to use both crop() and
mask() together in most cases. This combination of
functions would (a) limit the raster’s extent to our area of interest
and then (b) replace all of the values outside of the area to NA.
srtm_cropped = crop(srtm, zion)
srtm_final = mask(srtm_cropped, zion)
plot(srtm_final)
Changing the settings of mask() yields different
results. Setting updatevalue = 0, for example, will set all
pixels outside the national park to 0. Setting
inverse = TRUE will mask everything inside the bounds of
the park (see ?mask for details)
srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
plot(srtm_inv_masked)
Raster extraction
Raster extraction is the process of identifying and returning the
values associated with a ‘target’ raster at specific locations, based on
a (typically vector) geographic object location. The results depend on
the type of selector used (points, lines or polygons) and arguments
passed to the extract() function, which we use to
demonstrate raster extraction. The reverse of raster extraction —
assigning raster cell values based on vector objects — is
rasterization.
The basic example is of extracting the value of a raster cell at
specific points. For this purpose, we will use zion_points,
which contain a sample of 30 locations within the Zion National
Park.
# access some available data called zion_points
data("zion_points", package = "spDataLarge")
# now plot them all
plot(srtm)
plot(zion, add=TRUE, col=NA)
plot(zion_points, add=TRUE, col='red')
The following command extracts elevation values from
srtm and creates a data frame with points IDs
(one value per vector’s row) and related srtm values for
each point.
elevation = extract(srtm, zion_points)
elevation## [1] 1802 2433 1886 1370 1452 1635 1380 2032 1830 1860 1440 2145 1942 1691 1776
## [16] 2198 1820 1349 1758 1424 2159 1809 1826 1550 1799 2102 2118 1372 1905 1574Now, we can add the resulting object to our zion_points
dataset with the cbind() function:
zion_points = cbind(zion_points, elevation) # insert the elevation data into the orginal shapefile
zion_points## Simple feature collection with 30 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2077 ymin: 37.16632 xmax: -112.8717 ymax: 37.43165
## Geodetic CRS: WGS 84
## First 10 features:
## elevation geometry
## 1 1802 POINT (-112.9159 37.20013)
## 2 2433 POINT (-113.0937 37.39263)
## 3 1886 POINT (-113.0246 37.33466)
## 4 1370 POINT (-112.9611 37.24326)
## 5 1452 POINT (-112.9898 37.20847)
## 6 1635 POINT (-112.8807 37.19319)
## 7 1380 POINT (-113.0505 37.24061)
## 8 2032 POINT (-113.0953 37.34965)
## 9 1830 POINT (-113.0362 37.31429)
## 10 1860 POINT (-113.2077 37.43165)Note that now we have the elevation values for each one of the points
associated with zion_points. For fun let’s use
mapview to make an interactive map of both. Again we will
use zcol to set the column used for symbology in the point
map.
if(!require("mapview")) install.packages("mapview")## Loading required package: mapviewmapview(srtm)+
mapview(zion_points, zcol='elevation')Let’s tie this lesson back to our work with OSM. Here we are going to
import the package osmdata to allow us to interact with the
Overpass Turbo API. Let’s try to build a query using a bounding box from
getbb(), where we search for Zion National Park, this
bounding box then get passed to Overpass Turbo and restaurants will be
added.
# import functions from osmdata library
if(!require("osmdata")) install.packages("osmdata")
# create a overpass turbo query
query = getbb(place_name="Zion National Park, Utah, USA") %>%
opq() %>%
add_osm_feature(key="amenity", value="restaurant")
# convert to a spatial data type sf
restaurant = osmdata_sf(q=query)Ok, bad request! Thats not too specific. Let’s try to see what happened. First let’s see if we were able to get a bounding box for Zion National Park.
getbb(place_name="Zion National Park, Utah, USA")Looks like that is the problem. You can try to see if you can adjust
the name to get a valid response. But its likely we will need to create
a bounding box manually. Luckily we can get the bounding box from the
srtm raster file using bbox(srtm), this then
needs to be converted to a vector using c() in order to
meet the requirement of the opq function. Look at
?opq, notice that we can pass a bounding box in the format
bbox = c(xmin, ymin, xmax, ymax).
Note: bounding boxes can also be found using tools like bboxfinder.com.
Now let’s try again, using our new bounding box.
query = opq(bbox = c(bbox(srtm))) %>%
add_osm_feature(key="amenity", value="restaurant")
restaurant = osmdata_sf(q=query)
restaurant## Object of class 'osmdata' with:
## $bbox : 37.1320834298579,-113.239583212784,37.5129167631658,-112.85208321281
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 32 points
## $osm_lines : NULL
## $osm_polygons : 'sf' Simple Features Collection with 4 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULLLet’s just focus on the points that we returned to us:
restaurant_points = restaurant$osm_points
restaurant_points## Simple feature collection with 32 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2018 ymin: 37.16199 xmax: -112.8757 ymax: 37.27104
## Geodetic CRS: WGS 84
## First 10 features:
## osm_id name addr.city addr.housenumber
## 880118222 880118222 <NA> <NA> <NA>
## 883329842 883329842 Zion Pizza & Noodle Co. <NA> <NA>
## 904095151 904095151 Pioneer restaurant <NA> <NA>
## 937546158 937546158 Parallel Eighty-Eight <NA> <NA>
## 1579995974 1579995974 Thai Sapa Springdale 198
## 2983983311 2983983311 Switchback Grille Steakhouse <NA> <NA>
## 2983983665 2983983665 Arkansas Al's Steakhouse <NA> <NA>
## 3641061615 3641061615 <NA> <NA> <NA>
## 3641061616 3641061616 <NA> <NA> <NA>
## 3641061617 3641061617 <NA> <NA> <NA>
## addr.postcode addr.state addr.street addr.unit amenity
## 880118222 <NA> <NA> <NA> <NA> <NA>
## 883329842 <NA> <NA> <NA> <NA> restaurant
## 904095151 <NA> <NA> <NA> <NA> restaurant
## 937546158 <NA> <NA> <NA> <NA> restaurant
## 1579995974 84767 UT Zion Park Boulevard <NA> restaurant
## 2983983311 <NA> <NA> <NA> <NA> restaurant
## 2983983665 <NA> <NA> <NA> <NA> restaurant
## 3641061615 <NA> <NA> <NA> <NA> <NA>
## 3641061616 <NA> <NA> <NA> <NA> <NA>
## 3641061617 <NA> <NA> <NA> <NA> <NA>
## cuisine level name.en name.fr
## 880118222 <NA> <NA> <NA> <NA>
## 883329842 pasta;pizza 1 <NA> <NA>
## 904095151 <NA> <NA> <NA> <NA>
## 937546158 <NA> <NA> <NA> <NA>
## 1579995974 thai <NA> <NA> <NA>
## 2983983311 fish;steak_house <NA> <NA> <NA>
## 2983983665 american <NA> <NA> <NA>
## 3641061615 <NA> <NA> <NA> <NA>
## 3641061616 <NA> <NA> <NA> <NA>
## 3641061617 <NA> <NA> <NA> <NA>
## opening_hours
## 880118222 <NA>
## 883329842 Mo-Su 16:00-21:00
## 904095151 <NA>
## 937546158 <NA>
## 1579995974 <NA>
## 2983983311 Apr-Oct 17:00-21:00; Dec-Feb Mo-We off; Dec-Feb Th-Su 17:00-20:00; Mar,Nov 17:00-20:00; Dec 24-25 off
## 2983983665 <NA>
## 3641061615 <NA>
## 3641061616 <NA>
## 3641061617 <NA>
## outdoor_seating phone takeaway
## 880118222 <NA> <NA> <NA>
## 883329842 yes +1-435-772-3815 yes
## 904095151 <NA> <NA> <NA>
## 937546158 <NA> <NA> <NA>
## 1579995974 <NA> +14357720510 <NA>
## 2983983311 <NA> +1-435-772-3700 <NA>
## 2983983665 <NA> <NA> <NA>
## 3641061615 <NA> <NA> <NA>
## 3641061616 <NA> <NA> <NA>
## 3641061617 <NA> <NA> <NA>
## website geometry
## 880118222 <NA> POINT (-113.0003 37.18914)
## 883329842 http://www.zionpizzanoodle.com/ POINT (-112.9984 37.18924)
## 904095151 <NA> POINT (-112.9976 37.18926)
## 937546158 <NA> POINT (-113.0062 37.17994)
## 1579995974 http://www.thaisapazion.com/ POINT (-112.9918 37.19809)
## 2983983311 http://switchbackgrille.com/grille/ POINT (-113.0026 37.18402)
## 2983983665 <NA> POINT (-113.014 37.16757)
## 3641061615 <NA> POINT (-112.8759 37.27104)
## 3641061616 <NA> POINT (-112.8757 37.2709)
## 3641061617 <NA> POINT (-112.8758 37.27078)That is too much data. Let’s select only those columns
we are interested in:
restaurant_points = subset(restaurant_points, select = c(name,phone ))
restaurant_points## Simple feature collection with 32 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2018 ymin: 37.16199 xmax: -112.8757 ymax: 37.27104
## Geodetic CRS: WGS 84
## First 10 features:
## name phone
## 880118222 <NA> <NA>
## 883329842 Zion Pizza & Noodle Co. +1-435-772-3815
## 904095151 Pioneer restaurant <NA>
## 937546158 Parallel Eighty-Eight <NA>
## 1579995974 Thai Sapa +14357720510
## 2983983311 Switchback Grille Steakhouse +1-435-772-3700
## 2983983665 Arkansas Al's Steakhouse <NA>
## 3641061615 <NA> <NA>
## 3641061616 <NA> <NA>
## 3641061617 <NA> <NA>
## geometry
## 880118222 POINT (-113.0003 37.18914)
## 883329842 POINT (-112.9984 37.18924)
## 904095151 POINT (-112.9976 37.18926)
## 937546158 POINT (-113.0062 37.17994)
## 1579995974 POINT (-112.9918 37.19809)
## 2983983311 POINT (-113.0026 37.18402)
## 2983983665 POINT (-113.014 37.16757)
## 3641061615 POINT (-112.8759 37.27104)
## 3641061616 POINT (-112.8757 37.2709)
## 3641061617 POINT (-112.8758 37.27078)Now let’s plot it out and assign the restaurant names as labels:
mapview(srtm)+
mapview(restaurant_points,
label = restaurant_points$name,
legend=F)Let’s extract the elevation of each restaurant, following the same code as earlier:
elevation = extract(srtm, restaurant_points)
elevation## [1] 1197 1194 1194 1173 1212 1184 1175 1984 1984 1984 1984 1153 1153 1162 1162
## [16] 1210 1202 1200 1200 1202 1195 1194 1197 1184 1197 1197 1200 1089 1318 1192
## [31] 1191 1196Now, we can add the resulting list of elevation values to our
restaurant_points dataset with the cbind()
function:
restaurant_points = cbind(restaurant_points, elevation) # insert the elevation data into the orginal shapefile
restaurant_points## Simple feature collection with 32 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2018 ymin: 37.16199 xmax: -112.8757 ymax: 37.27104
## Geodetic CRS: WGS 84
## First 10 features:
## name phone elevation
## 880118222 <NA> <NA> 1197
## 883329842 Zion Pizza & Noodle Co. +1-435-772-3815 1194
## 904095151 Pioneer restaurant <NA> 1194
## 937546158 Parallel Eighty-Eight <NA> 1173
## 1579995974 Thai Sapa +14357720510 1212
## 2983983311 Switchback Grille Steakhouse +1-435-772-3700 1184
## 2983983665 Arkansas Al's Steakhouse <NA> 1175
## 3641061615 <NA> <NA> 1984
## 3641061616 <NA> <NA> 1984
## 3641061617 <NA> <NA> 1984
## geometry
## 880118222 POINT (-113.0003 37.18914)
## 883329842 POINT (-112.9984 37.18924)
## 904095151 POINT (-112.9976 37.18926)
## 937546158 POINT (-113.0062 37.17994)
## 1579995974 POINT (-112.9918 37.19809)
## 2983983311 POINT (-113.0026 37.18402)
## 2983983665 POINT (-113.014 37.16757)
## 3641061615 POINT (-112.8759 37.27104)
## 3641061616 POINT (-112.8757 37.2709)
## 3641061617 POINT (-112.8758 37.27078)Extract Values Along a Line
Raster extraction also works with line selectors. Then, it extracts
one value for each raster cell touched by a line. In this case, the best
approach is to split the line into many points and then extract the
values for these points. To demonstrate this, the code below creates
zion_transect, a straight line going from northwest to
southeast of the Zion National Park
Let’s create a line based on two points
c(-113.2, -112.9) and c(37.45, 37.2) which are
xy pairs in lat lon. Note that we are assigning a projection using an
EPSG code of 4326. EPSG codes are unique identifying codes
that are available for all common projections. To look up EPSG codes
please go to https://epsg.io.
zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
st_linestring() %>%
st_sfc(crs = 'EPSG:4326')%>%
st_sf()
zion_transect## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## geometry
## 1 LINESTRING (-113.2 37.45, -...Now let’s see if its showing up in the right place
mapview(srtm)+
mapview(zion_transect,
color='white')The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an ‘elevation profile’ of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.
Because our line is made up of only two points and
extract works best with points, we need to break our longer
line into lots of points.
The first step is to add a unique id for each transect. Next, with
the st_segmentize() function we can add points along our
line(s) with a provided density (dfMaxLength) and convert
them into points with st_cast().
Let’s look at st_segmentize help to see what units
dfMaxLength is in:
?st_segmentizeNow lets create our collection of points:
zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")
zion_transect## Simple feature collection with 257 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## First 10 features:
## id geometry
## 1 1 POINT (-113.2 37.45)
## 1.1 1 POINT (-113.1988 37.44902)
## 1.2 1 POINT (-113.1976 37.44805)
## 1.3 1 POINT (-113.1965 37.44707)
## 1.4 1 POINT (-113.1953 37.4461)
## 1.5 1 POINT (-113.1941 37.44512)
## 1.6 1 POINT (-113.1929 37.44415)
## 1.7 1 POINT (-113.1918 37.44317)
## 1.8 1 POINT (-113.1906 37.4422)
## 1.9 1 POINT (-113.1894 37.44122)Now, we have a large set of points, and we want to derive a distance between the first point in our transects and each of the subsequent points. In this case, we only have one transect, but the code, in principle, should work on any number of transects:
zion_transect = zion_transect %>%
group_by(id) %>%
mutate(dist = st_distance(geometry)[, 1])
zion_transect## Simple feature collection with 257 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## # A tibble: 257 × 3
## # Groups: id [1]
## id geometry dist
## * <int> <POINT [°]> [m]
## 1 1 (-113.2 37.45) 0
## 2 1 (-113.1988 37.44902) 150.
## 3 1 (-113.1976 37.44805) 300.
## 4 1 (-113.1965 37.44707) 450.
## 5 1 (-113.1953 37.4461) 600.
## 6 1 (-113.1941 37.44512) 750.
## 7 1 (-113.1929 37.44415) 901.
## 8 1 (-113.1918 37.44317) 1051.
## 9 1 (-113.1906 37.4422) 1201.
## 10 1 (-113.1894 37.44122) 1351.
## # … with 247 more rowsFinally, we can extract elevation values for each point in our transects and combine this information with our main object.
zion_elev = extract(srtm, zion_transect)
zion_transect = cbind(zion_transect, zion_elev)
zion_transect## Simple feature collection with 257 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## First 10 features:
## id dist zion_elev geometry
## 1 1 0.0000 [m] 2001 POINT (-113.2 37.45)
## 2 1 150.0961 [m] 2033 POINT (-113.1988 37.44902)
## 3 1 300.1922 [m] 2020 POINT (-113.1976 37.44805)
## 4 1 450.2883 [m] 1989 POINT (-113.1965 37.44707)
## 5 1 600.3844 [m] 1902 POINT (-113.1953 37.4461)
## 6 1 750.4805 [m] 1869 POINT (-113.1941 37.44512)
## 7 1 900.5766 [m] 1827 POINT (-113.1929 37.44415)
## 8 1 1050.6728 [m] 1767 POINT (-113.1918 37.44317)
## 9 1 1200.7689 [m] 1772 POINT (-113.1906 37.4422)
## 10 1 1350.8650 [m] 1722 POINT (-113.1894 37.44122)We can now see that we have both the distance from the first point in our line, and the elevation of each point along that line.
Let’s plot the transect and see what it looks like:
plot(x=zion_transect$dist,
y = zion_transect$zion_elev,
type='l',
col='darkred',
main='Elevation Profile',
xlab = 'Distance',
ylab = 'Elevation (m)')
Step 1 - Create a query to pull all the ‘primary’ roads in Zion
National Park and save them in a variable called roads.
Make sure to look at the tagging convention on OSM’s map
features directory. Also you need to set the bbox from
the query, you will want to do it based on the srtm raster
as we did earlier.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# # hint
# query = opq(bbox = ?????) %>%
# add_osm_feature(key="?????", value="?????")
#
# roads = osmdata_sf(q=query)Step 2: osmdata_sf returns points
(roads$osm_points), lines (roads$osm_lines),
and polygons(roads$osm_polygons). Create a new variable
called primary_roads that holds the road line features.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# print(primary_roads)Step 3: Use the filter() function from
dpylr to isolate only roads with the name
Zion Park Boulevard, and store the output in a variable
called park_road:
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# print(park_road)Let’s take a look at the data:
park_roadNotice that there are multiple segments making up this one road. Roads are typically broken into segments to allow for connectivity between other roads, and to maintain data like speed limits - which likely vary along the road.
Step 4: We need to ‘dissolve’ all of these road features into a
single line representing ‘Zion Park Boulevard’. In order to do this we
are going to tell the computer to create a new geometry
column using summarize, based on the st_union
of the multiple geometry rows. This is one that we haven’t talked about
yet, so I will just do it for you.
park_road = park_road %>%
summarize(geometry = st_union(geometry))
plot(park_road)
print(park_road)## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -113.0368 ymin: 37.16021 xmax: -112.9902 ymax: 37.19982
## Geodetic CRS: WGS 84
## geometry
## 1 LINESTRING (-113.0368 37.16...Notice there is now only one row of data and our multiple road segments are now the union.
Step 5: Now we can create our column called id, then use
st_segmentize to break the road into multiple points at a
fixed distance (250 or 500 works well), then we need to cast that point
object back into an sf point object using
st_cast.
See if you can find the correct chunk of code
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# park_road$id = 1:nrow(???)
# park_road = st_segmentize(park_road, dfMaxLength = ????)
# park_road = st_cast(park_road, "?????")
# park_roadStep 6: Let’s calculate the distances from the first point to all following points:
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# park_road = park_road %>%
# group_by(?????) %>%
# mutate(dist = st_distance(?????)[, 1])
# park_roadStep 7: Now we can use the extract function to grab the
values from the raster srtm at each one of the points
making up park_road, then we need to cbind
those values back into the park_road_elev dataset.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# park_road_elev = extract(????, park_road)
# park_road_elev = cbind(????, ????)
# park_road_elevStep 8: Plot the park_road_elev$dist on the x axis and
park_road_elev$park_road_elev on the y axis, and set the
type to 'line'.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# plot(x = park_road_elev$?????,
# y = park_road_elev$?????,
# type='line')